clc;
clear;
close all;

%%
%生成发生数据

M=4;
k=log2(M);
f=imread("coolLearn.PNG");
data=reshape(f,[],1);
bData=dec2bin(data,8);
bDate1=reshape(bData,[],k);
decData=bin2dec(bDate1);

%序列补位
N=length(decData);
zero_num=mod(5-mod(N,5),5);
decData_0=[decData;zeros(zero_num,1)];

%%
%QPSK调制
f_psk=pskmod(decData_0,M,pi/M,"gray");
X=reshape(f_psk,5,[]);      %将发送数据并行分成5个流

%设置信道参数
N_TX=5;
N_RX=5;
H=(randn(N_RX,N_TX)+1j*randn(N_RX,N_TX))/sqrt(2);  %瑞利衰落信道下的信道矩阵
SNRdB=0:2:10;   %信噪比

%初始化检测参数
berZF=zeros(length(SNRdB),1);%误码率
berSIC=zeros(length(SNRdB),1);
berML=zeros(length(SNRdB),1);
berMMSE=zeros(length(SNRdB),1);

timeZF=zeros(length(SNRdB),1);%时间
timeSIC=zeros(length(SNRdB),1);
timeML=zeros(length(SNRdB),1);
timeMMSE=zeros(length(SNRdB),1);

%%
% 天线接收
for i=1:length(SNRdB)
    %过AWGN信道
    sigma = sqrt(10^(-SNRdB(i)/10)); % 信噪比计算公式
    n = sigma*(randn(N_RX,length(X))+1i*randn(N_RX,length(X)))/sqrt(2);   
    % 生成接收信号
    Y = H*X + n;

    %%
    %线性检测——ZF符号检测 

    tstart=tic; %计算运行时间

    tZF=pinv(H.'*H)*H.'*Y;
    tZF1=reshape(tZF,[],1);
    sZF = pskdemod(tZF1,M,pi/M,'gray');
    toc(tstart);
    timeZF(i)= toc(tstart);%记录运行时间

    [~,berZF(i)]=biterr(decData_0,sZF) %计算误码率

    %恢复数据
    sZF=sZF(1:end-zero_num);
    ZFb1=dec2bin(sZF,2);
    ZFb2=reshape(ZFb1,[],8);
    ZFb3=bin2dec(ZFb2);
    ZFb4=reshape(ZFb3,size(f))/255;

    %画猫图
    figure(1)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(ZFb4));
    title(['ZF SNR=',num2str(SNRdB(i)),'dB']);


%%
% SIC检测
    tstart=tic; %计算运行时间
    squares_sum=sum(abs(H).^2,1);
    [~,sort_idx] =sort(squares_sum);
    H_sorted=H(:,sort_idx);
    Y2=Y;
    tSIC = zeros(size(X));
    for sic_i =1 :N_TX
         X_temp = pinv(H_sorted.'*H_sorted)*H_sorted.'*Y2;
         tSIC(end-sic_i+1,:) = pskdemod(X_temp(end,:),4,pi/M,'gray');
         Y2=Y2 - H_sorted(:,end)*pskmod(tSIC(end-sic_i+1,:),4,pi/M,'gray');
         H_sorted = H_sorted(:,1:end-1); %删除已估计的列
    end

     % 恢复原始矩阵行数顺序
    [~, original_idx] = sort(sort_idx);
    tSIC = tSIC(original_idx, :);
    toc(tstart);
    timeSIC(i) = toc(tstart); % 记录运行时间
    sSIC = reshape(tSIC, [], 1);
 
    [~,berSIC(i)] = biterr(decData_0 , sSIC)%计算SIC误码率
 
    %恢复数据
    sSIC = sSIC(1: end-zero_num);
    SICb1 = dec2bin(sSIC,2);
    SICb2 = reshape(SICb1,[],8);
    SICd3 = bin2dec(SICb2);
    SICd4 = reshape(SICd3,size(f))/255;
    
    %画猫图
    figure(2)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(SICd4));
    title(['SIC SNR=',num2str(SNRdB(i)),'dB']);

%% 
%ML检测 

   tstart = tic; % 计算运行时间
    permutation_list = (dec2base(0:M^N_TX-1, 4) - '0').'; 
    permutation_list_psk = pskmod(permutation_list,4, pi/M,"gray");  
    num_permutations = size(permutation_list,2);       
    distance = zeros(1,num_permutations);      
    tML = zeros(size(X));                      
    for ml_i = 1:length(Y)
        ml_xk_matrix = H*permutation_list_psk;
        distance = sum(abs(Y(:,ml_i) - ml_xk_matrix).^2);
        [~, ml_index] = min(distance);
        tML(:,ml_i) = permutation_list(:,ml_index);
    end

    toc(tstart);
    timeML(i) = toc(tstart); % 记录运行时间
    sML = reshape(tML, [], 1);             
    
    [~,berML(i)] = biterr(decData_0 , sML)  %计算ML误码率
 
    %恢复数据
    sML = sML(1: end-zero_num);           
    MLb1 = dec2bin(sML,2);
    MLb2 = reshape(MLb1,[],8);
    MLd1 = bin2dec(MLb2);
    MLd2 = reshape(MLd1,size(f))/255;
    
    %画猫图
    figure(3)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(MLd2));
    title(['ML SNR=',num2str(SNRdB(i)),'dB']);
%%
% MMSE
    tstart = tic; % 计算运行时间

    gamma = 10 ^ (- SNRdB(i) / 10); 
    LAMBDA = gamma * eye(N_RX); 
    wMMSE = H' * pinv(H * H' + LAMBDA);
    tMMSE = wMMSE * Y;
    tMMSE1 = reshape(tMMSE, [], 1);
    sMMSE = pskdemod(tMMSE1, M,  pi/M,'gray');
    toc(tstart);
    timeMMSE(i) = toc(tstart);
 
    [~,berMMSE(i)] = biterr(decData_0 , sMMSE);%计算误码率
   

    sMMSE = sMMSE(1: end-zero_num);
    MMSEb1 = dec2bin(sMMSE,2);
    MMSEb2 = reshape(MMSEb1,[],8);
    MMSEd1 = bin2dec(MMSEb2);
    MMSEd2 = reshape(MMSEd1,size(f))/255; 
    
    %画猫图
    figure(4)
    subplot(length(SNRdB),1,i);
    imshow(im2uint8(MMSEd2));
    title(['MMSE SNR=',num2str(SNRdB(i)),'dB']);
end

%绘制结果
figure(5)
semilogy(SNRdB, berZF, 'b-o',...
    SNRdB, berMMSE, 'r-s',...
    SNRdB, berSIC, 'g-d',...
    SNRdB, berML, 'm-p', 'LineWidth', 2);
grid on;
xlabel('SNR(dB)');
ylabel('BER');
legend('ZF', 'MMSE', 'SIC', 'ML');

figure(6);
semilogy(SNRdB, timeZF, 'b-o',...
    SNRdB, timeMMSE, 'r-s',...
    SNRdB, timeSIC, 'g-d',...
    SNRdB, timeML, 'm-p','LineWidth', 2);
grid on;
xlabel('SNR(dB)');
ylabel('Time (s)');
legend('ZF', 'MMSE', 'SIC', 'ML');
